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ABSTRACT 

We apply the iterative Expectation-Maximization algorithm (EM) to estimate the 
power spectrum of the CMB from multifrequency microwave maps. In addition, we are 
also able to provide a reconstruction of the CMB map. By assuming that the combined 
emission of the foregrounds plus the instrumental noise is Gaussian distributed in 
Fourier space, we have simplified the EM procedure finding an analytical expression 
for the maximization step. By using the simplified expression the CPU time can be 
greatly reduced. We test the stability of our power spectrum estimator with realistic 
simulations of Planck data, including point sources and allowing for spatial variation of 
the frequency dependence of the Galactic emissions. Without prior information about 
any of the components, our new estimator can recover the CMB power spectrum up 
to scales I « 1500 with less than 10 % error. This result is significantly improved if 
the brightest point sources are removed before applying our estimator. In this way, 
the CMB power spectrum can be recovered up to / « 1700 with 10 % error and up to 
I « 2100 with 50 % error. This result is very close to the one that would be obtained 
in the ideal case of only CMB plus white noise, for which all our assumptions are 
satisfied. Moreover, the EM algorithm also provides an straightforward mechanism to 
reconstruct the CMB map. The recovered cosmological signal shows a high degree of 
correlation (r = 0.98) with the input map and low residuals. 

Key words: cosmic microwave background, niethods:statistical 



1 INTRODUCTION 

Undergoing CMB experiments like BOOMERanG (Netter- 
field et al. 2002, Rulh et al. 2002), MAXIMA (Hanany et 
al. 2000), DASI (Halverson et al. 2002), VSA (Grainge et 
al. 2002), CBI (Mason et al. 2002), ACBAR (Kuo et al. 
2002), Archeops (Benoit et al. 2002) and MAP as well as 
future ones (Planck) will reveal with unprecedented quality 
the primordial matter fluctuations responsible for the CMB. 
Very recently the first detection of the E-mode polarization 
in the CMB has been claimed (DASI, Kovac et al. 2002), 
which independently supports the structure formation mod- 
els via gravitational instability. Thanks to these experimen- 
tal results the fundamental cosmological parameters can be 
determined with good accuracy. 

These experiments will measure also the emission com- 
ing from our own Galaxy at mm frequencies as well as the 
emission due to other galaxies in the same wavebands. On 
the other hand, galaxy clusters will distort the CMB radi- 



ation with an intensity which is proportional to the total 
mass of the cluster times its temperature. All these compo- 
nents will generate a confusion limit which, together with 
the intrinsic detector noise, will make very diflicult to dis- 
entangle which percentage of the observed emission is due 
to the CMB or to the other components. Several component 
separation methods have been proposed so far in the litera- 
ture: multifrequency Wiener filter (MWF, Tegmark and Efs- 
tathiou 1996, Bouchet and Gispert 1999), maximum entropy 
methods (MEM, Hobson et al. 1998, Vielva et al. 2001b; ex- 
tended to the sphere in Stolyarov et al. 2002, Barreiro et al. 
2003), independent component analysis (ICA, Baccigalupi 
et al. 2000, Maino et al. 2002), blind Bayes (Snoussi et al. 
2001). 

Although there are some similarities and differences between 
the previous methods, all of them share a common aspect: 
they try to separate all the components simultaneously. For 
this purpose, these methods usually need to assume some a 
priori information about the components. Thus, if the differ- 
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ent components have different frequency dependencies, then, 
by examining the data at different frequencies it is possible 
to distinguish (at a certain level) the different components. 
Or, if we know the correlation function of the components 
(or their power spectrum) and each one is significantly dif- 
ferent from each other, then, it is also possible to distinguish 
(again at a certain level) the components since their modes 
will behave differently in the Fourier space. If one knows 
both, the frequency dependence and the power spectrum, 
then the component separation improves dramatically since 
now it is possible to combine the information in the different 
channels by correlating the Fourier modes with a correlation 
given essentially by the frequency dependence of each com- 
ponent and its power spectrum. This approach has proved 
to be very useful if the frequency dependence and power 
spectrum of the components are known. 
It is interesting to explore what kind of information can be 
obtained when no prior information is assumed about none 
of the components. In some cases wc simply do not know the 
priors (as it happens in the case of the free-free emission or 
the spinning dust). On the other hand, if we assume some- 
thing about the components which is not accurate but an 
approach, then we are intrinsically introducing some bias in 
our final result. If, after the component separation, we end 
up with a biased estimate on, at least, one of the compo- 
nents, this will have an effect on some (if not all) of the 
other components which must absorb that bias in order to 
obey the constraint that the sum of all the components must 
be equal to the original data. This very last point is one of 
the risks of the simultaneous component separation meth- 
ods. 

Another assumption usually made is that the components 
arc, one from each other, statistically independent. This as- 
sumption is not true for several of the components (like the 
Galactic ones). This assumption can be a source of system- 
atic errors. Other typical assumption is that the probabil- 
ity distribution function (pdf) of the individual components 
is a Gaussian. This is false for components like the point 
sources, the Sunyaev-Zeldovich (SZ) and the Galactic emis- 
sions where the true pdf has a bell-shape with a long tail 
towards positive values (negative for the SZ for frequencies 
below 217 GHz). This assumtion will be somewhat relaxed 
in this work by assuming that the combined foreground 
emission plus instrumental noise is Gaussian distributed in 
Fourier space, instead of assuming it for each individual com- 
ponent. Wc will study the effect of this assumption at dif- 
ferent scales. At those scales where the assumption is less 
justified solutions will be proposed to reduce the effect. 

In a recent paper (Delabrouille et al. 2002) the authors 
have used the Expectation-Maximization (EM) algorithm 
to minimize a spectral matching criterion. They consider 
the problem of component separation with four components, 
CMB, SZ effect, dust and instrumental noise. In that work, 
the authors have shown that this is an interesting alternative 
to estimate jointly the frequency dependence and spatial 
power spectra of the components. This method allows to 
introduce certain degrees of freedom or free parameters in 
the problem. These free parameters can be determined after 
iterating an expectation-maximization process. 

In the work presented here, we will explore this direc- 
tion but in a much more simplified manner focusing only 
on the estimation of the power spectrum and the map re- 



construction of the CMB but including more components 
than in the previous work and doing no assumptions at all 
about the power spectrum or frequency dependence of any 
of the components. In addition to the diffuse foregrounds, 
we also will study the effect of point sources which were not 
considered in the previous work. 

The approach presented in this paper can be extended 
to determine jointly the CMB and the SZ components. The 
physics of these two emissions is well known, in particu- 
lar their frequency dependence. Therefore, no prior assump- 
tions about their frequency behaviour is required. In addi- 
tion, these two components -together with the point sources 
one- arc the most important signals from the cosmological 
point of view. Several works have already been presented 
to determine the SZ (Herranz et al. 2001b, c; Diego et al. 
2002) and the point sources emissions (Tegmark & Oliveira- 
Costa 1998; Cayon et al. 2000; Sanz et al. 2001; Vielva et al. 
2001, 2003; Chiang et al. 2002; Hobson & McLachlan 2002) 
from microwave images. Our aim for the future will be to de- 
velop a method based on the EM algorithm combined with a 
point source detection technique, to recover simultaneously 
the three cosmological emissions. 

The paper is organized as follows. In Section 2, wc sum- 
marize the key ideas behind the EM algorithm. In Section 3 
we apply the EM algorithm to the problem of determining 
an estimator of the CMB power spectrum. The main results 
are shown in Section 4. These results are compared with the 
ones which would be obtained in the ideal case when only 
CMB and white noise are considered. Finally, the conclu- 
sions are given in Section 5. 



2 THE EXPECTATION-MAXIMIZATION 
ALGORITHM 

In this section wc will present an alternative method which 
will be useful to get a robust estimate of the power spec- 
trum of the CMB. We will also present a single iterative ex- 
pression for the new estimator of the CMB power spectrum. 
This estimator can be used directly with multifrequency mi- 
crowave data to give a fast, accurate and robust estimate of 
the power spectrum of the CMB. 

EM is an algorithm useful when there is a many-to-one map- 
ping. That is, several variables are combined together to give 
one observed quantity. This is exactly our problem where the 
complete data are just the different components (Galactics, 
extragalactics, CMB and noise) and the observed quantity is 
just the projection of all of them in one of the channels (data 
on the receiver). EM allows to introduce degrees of freedom 
(or free parameters) in the pdf of the complete data and 
then look for the best parameters by maximising the ex- 
pected value of that pdf given the observations. 
The advantage of this algorithm in our problem is obvious. 
Since part of the information about the components is un- 
known, we can parametrise this unknown information as free 
parameters in the pdf of the complete data. Then, the max- 
imization process will give us the best set of free parameters 
given the data. 

The EM algorithm (Dempster et al. 1977) provides an 
iterative procedure for computing max;imum likelihood esti- 
mates (MLEs) in situations where the observed data vector, 
d, is viewed as being incomplete, d is an observable func- 



© 0000 RAS, MNRAS 000, 1-10 



CMB power spectrum with EM 3 



tion of the so-called complete data x (see e.g. McLachlan 
and Krishnan 1997). Let's denote by fi(d[p) the pdf's of 
the incomplete data d and by f(x|p) the pdf of x. The vec- 
tor p = (pi, pn) will contain the N unknown parame- 
ters. The complete-data log likelihood function that could 
be formed for p if x were fully observable is given by 



logL(p) = logf(x|p) 



(1) 



Instead of observing x, we observe d, with many-to-one map- 
ping from x to d {i.e d = d(x)). It follows that 



fi(d|p) 



f(xjp)dx 



(2) 



where Xj is the sub-space of x defined by d = d(x). 

The EM algorithm approaches the problem of solving the 

incomplete-data likelihood, log£(p) = logfi(d|p), indirectly 
by proceeding iteratively in terms of the complete-data log 
likelihood function, logLijp). As it is unobservable, it is re- 
placed by its conditional expectation given d, using the cur- 
rent fit for p. More specifically, the EM algorithm consists 
of two steps: Expectation (E-step) and Maximization (M- 
step). On the j-|-l iteration, these steps are as follows: 

• E-step: Calculate the quantity Q(p|p)^), defined as 
Q(p|p-') = i?{logL(p)|d,pJ}. 

• M-step: Choose p-""*"""" to be the value that maximises 
Q(p|pJ); that is, QCpS+^lpJ) = max[Q(p|pJ)]. 

That is, in the first step we compute the expected value 
of the log of the complete pdf given the data and an estimate 
of the parameters (f)^). In the second step we look for the 
values of the parameters which maximise the previous ex- 
pected value. At this step, we have to maximise with respect 
to each one of the parameters (p-' (i)). This can be done by 
just setting the first derivative of Q(p;p'') with respect to 
that parameter to and solving for 

This is an iterative process which can be started with an 
arbitrary choice for the parameters in the first step, p°. The 
EM algorithm assure us that at every step the likelihood 
£(p) is increased. The method will converge to the optimal 
set of parameters, p, in a number of steps which will depend 
on the nature of the problem and on the initial election for 
the first iteration, p°. 

In this paper, we will consider a simple case where wc 
allow the covariance in Fourier space of the CMB signal (or 
power spectrum) to be a free parameter. 
Wc will apply the EM algorithm in order to determine that 
free parameter (the power spectrum of the signal). This case 
can be extended to include more free parameters in the anal- 
ysis as it is shown in Delabrouille et al. (2002). 



3 EM APPLIED TO CMB EXPERIMENTS 

If we know the frequency dependence of one of the compo- 
nents, then we can express the signal at a given frequency 

as: 



S(n, v) = A(n, v) (g) (f{v)s{n)j 
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(3) 



where s(n) is the spatial pattern of the signal we want to 
estimate, f{i') is the frequency dependence of the signal (in- 
cluding the band-with), and A(n, i/)® accounts for the con- 
volution with the antenna beam of the experiment and its 

frequency response. 

The data at different frequencies, d(j/), can be expressed as 
a sum of the signal plus some other contributions. 



d(n,z/) = S(n,i.)-Fe(n,z/). 



(4) 



The residual ^(n, i^) includes all the other components, i.e. 
Galactic and extragalactic components convolved with the 
corresponding beam plus the corresponding instrumental 
noise. Due to the antenna and frequency convolution in 
A(n, I')®, it is easier to work in Fourier space where this con- 
volution is just a product of the Fourier modes which are un- 
correlated provided the field is homogeneous and isotropic. 
Therefore, the previous equation should be rather expressed 
in Fourier space as: 



(5) 



If the experiment has m diff'crcnt channels, then we have an 
observation for each channel. That is, we have the vector of 
observations d = {di,d2, ■■■,dm) for each Fourier mode k. 
And we can write: 



d = Rs + ^ 



(6) 



where, R is the response vector containing m elements. Each 
element of R is just the product of the antenna and the 
frequency response of the instrument times the frequency 
dependence of the signal. 



(7) 



In the case in which the signal is the CMB then fu = l V ^. 
From eqn. 6 we can express the residual as 



^ = d - Rs. 



(8) 



The basic starting point when applying EM is to define 
the pdf of the complete data. In our case the complete data 
are all the components plus the multifrequency observations 
d. We will focus on the problem in which the complete data 
can be divided in just two elements. The two elements are 
the CMB signal we want to estimate s and the observed data 
d (or equivatently the residual 5, noise plus rest of compo- 
nents; see below). By doing this, the nature of our problem 
reduces drastically since we only have to make assumptions 
about two elements. 

Furthermore, the residual and the signal can be considered 
as independent. In terms of the probability, the pdf of the 
complete data is: 



(9) 



f(x|p) = P(s,dlp) = P(s,^1p) = P(s1pa)P(^|pb) 

where p = IJpaPb with pA and pB unknown parameter 

sets for the pdf's of s and ^. The second equality follows 
from the transformation (s, d) — > (s,^) whose Jacobian is 
equal to one (see equation 5). It is important to note that, 
in oppossition to other methods, our assumption about in- 
dependence of the elements (s and ^) is well justified. No 
relation is expected between the CMB and the other fore- 
grounds (except maybe the SZ effect which could be weakly 
related to the CMB through the ISW effect). 
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Once we have established that the two elements are in- 
dependent, we need a clue about the specific form of the 
individual probabilities, P(s|pa) and P(^|pb). If the CMB 
is close to a Gaussian variable (and we know it is) , then the 
pdf of the CMB in Fourier space can be approached by: 



P(s|pa) = P(s) oc P{ky^exp{- 



P{k)' 



(10) 



where P{k) is the power spectrum of the CMB: the param- 
eter to be estimated. 

The pdf of the residual (^) can be modelled as 



P(e|pb) = P(0 oc exp{-iC-'i^) 



(11) 



where C^^ is the inverse of the correlation matrix of the 
residual. C is an n x n matrix (n = number of channels). 
This correlation matrix is not diagonal if there are some 
correlations in the residuals between the different frequency 
channels. The elements of C are the power (diagonal) and 
cross-power (off- diagonal) spectra of the residuals at each 
channel which can be given as a function of the P{k) (see 
below) . 

Wo have assumed in equation (11) that the residuals are 
Gaussian distributed. This assumption is not as strong as 
assuming that each component is Gaussian distributed. The 
sum of several non-Gaussian distributions tends to a more 
Gaussian one. In the limit of a sum of infinite independent 
components the central limit theorem assures that the re- 
sulting distribution is Gaussian. In any case Gaussianity can 
be a good approach at scales where the instrumental noise 
dominates the residuals. The most critical situations for the 
validity of this assumption appear at large and small scales 
since, at least for some channels, they arc dominated by the 
Galactic components and the compact sources, respectively. 
Solutions for these problems will be proposed in the next 
sections. 

In the previous relations for the individual pdf 's appear 
the power spectra of the signal and the residual. In a real 
situation these quantities are unknown. This is the reason 
why we apply EM. We can introduce the unknown CMB 
power spectrum, P{k), as the uncertainty in our problem. 
In order to apply EM and estimate this free parameter we 
have to calculate the expected value Q(p|p') which is just 
an integral: 



-j 



Q(p|pJ)= / /offP(s,d|p)P(s|d,p-')ds. 



(12) 



The only thing we need to know to calculate Q(p|i>^) are the 
probabilities P(s|d) and P(s,d). By Bayes we know that: 



P(s|d) = 



P(g,d) 
P(d) 



(13) 



where P(d) can be obtained just by marginalising P(s, d) 
over s: 



(14) 



P(d) = j V{s,d)ds 

and P(s, d) can be obtained through the chain rule: 



P(s,d) = P(s,C)J 



s, d 



(15) 



By taking equation (5), one can prove that the Jacobian in 
the previous expression is equal to one, as it was commented 
above, so we have: 



P(s,d)=P(6-,0 

where (see equation 9), 

2 

P(s,0 oc P(A;)"^ea;p( 



P{k) 



)exp{-iC-'e) 



(16) 



(17) 



where the term, P{k), accounts for our free parameter in 
P(s|pa) (see equation 10) and the second term was taken 
from equation (11). We remark that the covariance matrix 
C also depends on the free parameter P{k), as it will be 
discussed below. 

Now it is easy to compute the terms, P(d), P(s|d), and 
finally Q{P{k)\P{ky): 



P(d) oc 



exp(^ - dC-^t + (dC-'Rl')VFi^ 



P(fe)|C|l/2Fi 

where we have expressed ^ in terms of d and s (see eqn, 

Hr'~iRt 2 

P(s|d) oc Fiexp { -Fi 



(18) 



(19) 



where, 

Fi = RC-^Rt + P{k)-\ 



(20) 



Finally the expected value of the log of the complete pdf 
(equation 12) is: 

Q{P{k)\Piky) = -log{P{k))-^log\C\-dC-'d^ (21) 
- Fi < \sf > -|-2dC"^R''' < s > -l-const. 



where 



< s >= 



Fi 



, ,2 1 A |dC-iRt|2 



(22) 
(23) 



Equation (22) is nothing more than multifrcqucncy Wiener 
filter computed with the CMB power spectrum at iteration 
j. This is an expected result since Gaussianity for the sig- 
nal and residuals have been assumed. Similarly, the signal 
variance at each k < \sf > is also computed with P{ky 
at iteration j. At this point it is important to note that in 
eq. (22) the dependence on P{k) appears explicitly in the 
first term and also in the fourth one through the function 
Fi. Moreover, since we do not know the residuals the P{k) 
dependence is also present in the other terms through the 
residual covariance matrix, 



(24) 



where P^'^{k) and P^'^{k) are the residual and data cross- 
power spectra after convolution with the corresponding an- 
tenna response whereas P{k) is the unconvolved CMB power 
spectrum to be determined. 

In the maximisation step, the value of P{ky'^^ that 
maximises Q{P{k)\P{k)-') should be found. Finding the 
maximum is not a trivial task given the complicated depen- 
dence on P{k). In order to facilitate the finding we consider 
the residual covariance matrix, C^^ , constant in the M-step 
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(but changing with P{k) in the E-step). With this simpli- 
fication, wc can easily find the analitical solution for the 
maximum by differentiating Q(P(k)\P{ky) with respect to 
P{k) and equating to zero. We find the not so surprising 
result: 

P(fey+i =< \sf > . (25) 

The previous equation is the main result of EM. < |sp > 
is given by equation (23) computed with the value of P{k) 
at iteration j. By iterating equation (25), which only de- 
pends on the data, the response vector and the CMB power 
spectrum obtained in the previous iteration, we can get an 
estimate of the power spectrum of the underlying signal. 
The final step after the EM estimation of the signal power 
spectrum at each fc, is to recover the CMB map by applying 
eqn. 22, the multifrequency Wiener filter derived within the 
EM framework. 

Finally, this method can be easily extended to deter- 
mine jointly the CMB and the SZ emissions, since they are 
almost uncorrelated and the frequency dependence of the 
SZ is also very well known. In addition, the method can 
be implemented with more realistic (non-Gaussian) pdfs to 
model the signal distributions as well as the residuals. We 
will study these extensions of the method in a future work. 



Frequency 
(GHz) 


FWHM 
(circniin) 


Pixel size 
(circmiii) 




857 


5.0 


1.5 


22211 in 


545 


5.0 


1.5 


489.51 


353 


5.0 


1.5 


47.95 


217 


5.5 


1.5 


15.78 


143 


8.0 


1.5 


10.66 


100 (HFI) 


10.7 


3.0 


6.07 


100 (Lfl) 


10.1) 


:-!.() 


1 i..32 


70 


14.0 


3.0 


16.81 


44 


23.0 


6.0 


6.79 


30 


33.0 


6.0 


8.80 



Table 1. Experimental constrains at the 10 Planck channels. The 
antenna FWHM is given in column 2 for the different frequencies 
(a Gaussian pattern is assumed in the HFI and LFI channels). 
Characteristic pixel sizes are shown in column 3. The fourth col- 
umn contains information about the instrumental noise level, in 
AT/T per pixel. 



4 RESULTS 

In this section, we present the results that have been ob- 
tained for the CMB power spectrum determination. Firstly, 
we briefly describe the simulated data that have been used 
in this work, accounting for the most important contribu- 
tions to the microwave sky. Then, we present our CMB 
power spectrum estimation. An improvement at small scales 
is achived by using a very useful tool for subtracting the 
brightest point sources: the Mexican Hat Wavelet. 

Finally, we applied our estimator to an ideal data set, 
where the Gaussian assumption for the residuals is satisfied. 
This exercise is highly interesting, since it shows the robust- 
ness of our estimator: the results achieved in this idealistic 
case are only slightly better than the ones obtained in the 
reslistic situation. 

4.1 Simulated data 

To show the power of our approach we will apply our esti- 
mator (eqn. 25) to realistic simulations. They consider the 
expected levels of the instrumental white noise and reso- 
lution characteristics of the 10 Planck mission channels at 
30, 44, 70, 100 (Low and High Frequency Instrument, LFI 
and HFI), 143, 217, 353, 545 ad 857 GHz (see Table 1). 
In addition, we use state of the art simulations of the dif- 
ferent galactic and extragalactic components, together with 
the pure CMB signal. 

Within the first ones, we have taken into account the 
synchrotron, free-free, spinning and thermal dust contribu- 
tions. The first one has been simulated using the all-sky pat- 
tern provided by Giardino et al. (2002) including both, tem- 
perature and spectral index templates. The free-free emis- 
sion is very poorly known. Recent experiments focusing on 
the H-a detection, like Southern H-a Sky Survey (SHASSA, 
Reynolds & Haifner, 2000) and the Wisconsin H-a Mapper 



project (WHAM, Gaustad et al. 2001), will produce all-sky 
surveys that could be used as free-free templates. However, 
since these data axe not available at the present time, we 
have used the correlation between dust and free-free emis- 
sions proposed by Bouchet et al. (1996) as spatial template, 
and a power law, Ii, oc , to describe the frequency 

dependence. We have also included the spinning dust emis- 
sion -proposed by Draine & Lazarian (1998) as a possible 
explanation for the anomalous Galactic emission found by 
CMB experiments hke COBE (Kogut 1999) or Saskatoon 
(Oliveira-Costa et al. 1997). Very recently, a tentative con- 
firmation of that emission has been claimed (Finkbeiner et 
al. 2002). The simulation of that emission has been carried 
out using the frequency dependence proposed by Draine & 
Lazarian (1998)* and using the thermal dust component as 
an spatial template, since both dust emissions are strongly 
correlated through the neutral hydrogen column density 
(Nh)- Finally, thermal dust emission has been also simu- 
lated, using the best model proposed by Finkbeiner et al. 
(1999) to fit the FIRAS, IRAS and DIRBE data, consisting 
of two grey-bodies with mean temperatures of 16.2 K and 
9.4 K and emissivities 2.70 and 1.67, respectively. 

The simulated extragalactic foregrounds are the ther- 
mal Sunyaev-Zel'dovich effect (SZ) and radio and infrared 
point sources. The SZ was performed following the Diego et 
al. (2001) model for a flat ACDM Universe with Qrn = 0.3 
and Qa = 0.7. The point sources correspond to the Toffolatti 
et al. (1998) model for the same Universe, including radio 
fiat-spectrum and infrared sources. We refer to that paper 
for details. 

Finally, the CMB signal was simulated using the Ci's 
provided by the CMBFAST code (Seljak & Zaldarriaga, 
1996) for the same cosmological model and assuming Gaus- 

* Data provided by the authors. 
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S&SWse-: ^j^j. 




Figure 1. Realistic simulations of the data, covering a 12.8° x 12.8° patch, expected from the 10 Planck channels. The RMS amplitude 

(at 353 GHz and in AT/T thermodinamic temperature units) is 4.26 X lO"'"' for the CMB, 5.16 X 10"^ for the PS emission, 3.93 X 10""^ 
for the SZ contribution, 5.82 X 10~^ for the thermal dust, 3.23 X lO"'^ for the free-free radiation, 1.46 x 10^*' for the synchrotron emission 
and it is almost negligible in the case of the spinning dust (see text). 



© 0000 RAS, MNRAS 000, 1-10 



CMB power spectrum with EM 7 




Figure 2. Recovered power spectrum by our estimator. The solid 
line is the true power spectrum and the dotted and dashed lines 
represent the recovered power spectrum obtained by our estima- 
tor directly from the data and the one recovered after previous 
subtraction of the brightest point sources, respectively. P(k) is in 
AT/T units. Notice that the multipole / is related to the wave 
number k through I Ri 28fc. 
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k 

Figure 4. Relative error in % of the recovered power spectrum by 
our estimator after the brightest point sources have been removed 
in all the channels (dash-dotted line) . For comparison in the solid 

line we have plotted the errors obtained by our estimator in the 
ideal situation when only CMB and instrumental white noise are 
present in the data (see text). 
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Figure 3. Relative error in % of the recovered power spectrum by 
our estimator (dash-dotted line) . For comparison we have plotted 
the errors obtained by our estimator in the ideal situation when 
only CMB and instrumental white noise are present in the data 
(see text). 



sian fluctuations. We would like to remark that synchrotron, 
thermal dust and the point source emissions have been sim- 
ulated using a frequency dependence which varies with the 
sky position. The simulated maps covering a patch of the 
sky of 12.8° X 12.8°, as seen by the 10 Planck channels, are 
shown in figure 1. 

4.2 A simple choice for the correlation matrix 

In order to estimate the power spectrum of the CMB 
through cqn. (25) wc still have to iterate a relatively com- 
plex expression containing a cocient of polynomials of order 
twice the number of channels and involving many terms. To 
further simplify the convergence of that equation we will 



only consider the diagonal terms of the covariance matrix, 

making zero the other non-diagonal elements. This assump- 
tion does not take into account the correlations among the 
residuals at different channels due to the Galactic compo- 
nents and compact sources. The former are reflected at the 
low fe-modes and their effect can be reduced by performing 
a simple Galactic subtraction as the one proposed by Diego 
et al. (2002). The later have their impact at intermediate 
(e.g. radio sources observed in the lowest resolution chan- 
nels) and large A;-modes (e.g. infrared sources). However, as 
it will be shown below, once the brightest point sources are 
subtracted from the maps, the contribution of the remanent 
point sources is low enough to make tfieir effect significantly 
smaller. 

In figure 2 we present our results for the 10 Planck chan- 
nels using the data obtained from the combination of the 

simulated components as described above (figure 1). The 
comparison of the recovered spectrum with the true one 
shows a good agreement until k values larger than the one 
corresponding to the best antenna channels and the results 
start to be dominated by the deconvolved noise. 

In order to better appreciate the differences, we show in 
figure 3 the relative error, {Pr{k) — Pt{k))/Pt{k) %, where 
Pr{k) is the estimator of equation 25. 

In both large and relatively small scale range (small 
and relatively large k) our estimator finds a reasonably good 
estimation of the true power spectrum. The error is small 
(less than 10 %) up to scales fc w 55 (which corresponds 
to I « 1500). The effect of only considering the diagonal 
terms of the residual covariance matrix is small (see also the 
comparison with the ideal case below). This fact proves the 
applicability of our assumption for the microwave data since 
our method is not very sensitive to the cross-correlations of 
the different channels. 

At very small scales {k > 45 or Z > 1300) the signal 
is below the noise level and the error bars grow quickly. 
Our estimator shows a clear bias toward higher values (we 
recover more power than the true one). This implies that 
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some of our assumptions are wrong at these small scales. The 
Gaussianity assumption is a good approximation in Fourier 
space at large A;-modes (because of the dilution in many 
modes of the clear non-Gaussian features appearing in the 
real space). However, the uncorrclation assumption among 
frequency channels is a source of bias at these scales, since 
the bright point sources are strongly correlated. We have 
checked this point by removing the brightest point sources 
(but leaving the galaxy clusters and the weak point sources 
in the residual). Wo have used the Mexican Hat Wavelet 
(MHW) technique proposed by Vielva et al. (2001a) for re- 
moving point sources above the fluxes determined by the so- 
called 50% error criterion used in that paper. These fluxes 
correspond to those detection limits above which the max- 
imum percentage of spurious detection is 5%, being a de- 
tection spurious if the error in the amplitude estimation is 
larger than 50%. After the subtraction, we applied our esti- 
mator to these new point-source-free maps. The final result 
is shown in figure 4. The effect of removing the brightest 
point sources is evident now. This result clearly shows how 
the bias can be corrected by removing the compact sources. 
The relative error does not change significantly after remov- 
ing the point sources but the estimator is able to go up to 
scales = 60 (i « 1700) with relative error < 10%. The 
error is smaller than 50% below = 75 (Z « 2100). 

4.3 Comparing with an ideal situation 

It is interesting to answer whether or not the above results 
arc close to the most optimistic case where all our assump- 
tions are right (Gaussianity and C diagonal). We will check 
this point in the simple (but unrealistic) case where our data 
consist only of CMB signal plus the corresponding instru- 
mental noise for each channel. In this simple case, our as- 
sumption of Gaussianity is right since both CMB and noise 
are Gaussian. The correlation matrix of the residual, C, will 
be exactly diagonal with the elements being the power spec- 
trum of the data minus the CMB power spectrum convolved 
with the antenna corresponding to each channel. 

In figures 3 and 4 we show the performance of the EM 
estimator in this ideal case (CMB plus noise) after combin- 
ing the 10 Planck channels. The result is surprisingly very 
similar to the one obtained in the realistic case where the 
brightest point sources have been subtracted (figure 4). This 
test shows how the EM estimator is extremely robust and 
can give results very close to the most optimal one. 

4.4 Map reconstruction 

Using the expression given in equation 22 and the estimated 
power spectrum Pr{k), we are able to recover the CMB map. 
As it was pointed out in Section 3, due to the Gaussian- 
ity assumptions, equation 22 is nothing more than multifre- 
quency Wiener filter (MWF). Hence, our map reconstruc- 
tion is a MWF focused on the CMB recovery. Let us re- 
mark that this approach differs from the traditional MWF 
already used in the microwave sky recovery (e.g., Tegmark 
& Efstathiou 1996, Bouchct & Gispcrt 1999 and Hobson et 
al 1998). These works were focused on the all components 
separation, whereas our approach is devoted to reconstruct 
just the CMB signal. Even more, whereas the traditional 



approaches require previous knowledge of the power spec- 
tra (not only the CMB one, but also the power spectra of 
the foregrounds), the CMB power spectrum is a result of 
the blind EM method. This and the non requirement of any 
prior knowledge about the frequency dependence of the com- 
ponents, are the strongest points of the blind methods. In 
figure 5 we present the CMB map reconstruction (left) to- 
gether with the input one (right). Both maps have been con- 
volved with the best Planck resolution (FWHM = 5'). The 
correlation between the input and recovered CMB maps is 
very good r = 0.98 (figure 6), whereas the slope of the best 
straightline fit is 0.97. The residual CMB is dominated by 
the Galactic dust emission, since it is the most important 
foreground at low fc-modes. Hence, by performing a simple 
dust subtraction as the one proposed by Diego et al. (2002) 
-which basically consists in using the 857 GHz channels as 
a dust template to be subtracted from the others channels- 
a significant improvement in the CMB map reconstruction 
is achieved, being the RMS of the residual map lower than 
15%. 



5 CONCLUSION 

We have presented an efficient estimator (eqn. 25) of the 
CMB power spectrum. Our estimator is based on the EM 
algorithm and does not make use of any prior information 
whatsoever about any of the components. This renders sat- 
isfactory results, since we are able to recover the power spec- 
trum up to Z « 1500 with less than 10% error. This limit can 
be improved if we remove the brightest compact sources, for 
instance applying the MHW technique poposed by Vielva 
et al. (2001a). After point source removal the estimator can 
recover the power spectrum up to Z w 1700 with less than 
10% error and up to Z « 2100 with less than 50% error. 

Our results arc very close to the optimal case when all 
our assumptions are satisfied (only CMB and instrumen- 
tal noise are included in the simulations), as can be seen 
in figure 4. There are several CMB power spectrum estima- 
tions in the literature, given for different component separa- 
tion techniques, like Bayesian methods (MEM: Hobson et al. 
1998 and Stolyarov et al. 2002; and MWF: Tegmark & Efs- 
tathiou 1996 and Bouchet & Gispert 1999) and blind source 
separation algorithms (ICA and FastICA: Baccigalupi et 
al. 2000 and Maino et al. 2002; and Multi-Detector Multi- 
Component analysis: Delabrouille et al. 2002). This work 
provides a CMB power spectrum that is comparable to or 
even better than the previous ones, being more robust than 
those, since previous knowledge about the underlying com- 
ponents is not requiered (contrary to the Bayesian meth- 
ods) and simulations used in previous works arc significantly 
more idealised than the ones used in this work (e.g., lack 
of spatial variation of the frequency dependence for sev- 
eral components; absence of some of the components, like 
point sources or some Galactic foregrounds). On the other 
hand, our simulations account for the previous limitations 
including all the major microwave emissions and allowing for 
spatial variation of the frequency dependence. Finally, our 
method is able to reconstruct the CMB map using the esti- 
mated power spectrum and the MWF given in equation 22. 
Let us remark that the MWF is obtained when the Gaus- 
sianity hypothesis is assumed, whereas the EM algorithm 
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Figure 6. Correlation between the input and recovered CMB 
maps. The best strightline fit is also plotted. 



provides a more general proeedure for component separa- 
tion. That hypothesis will be explored in a future work by 
including more realistic distributions for the different com- 
ponents. In addition, we are working on a jointly deter- 
mination of the cosmological signals (CMB, SZ and point 
sources) using the EM algorithm and a point source de- 
tection tool. Finally, an extension of this algorithm to the 
sphere is straightforward and will be studied soon. 
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